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Numerical simulation of boundary layers 


Part 1. Weak formulation and numerical method 


By PHILIPPE R. SPALART 


NASA Ames Research Center, Moffett Field, California 94035 


A numerical method designed to solve the time-dependent, three-dimensional, incom- 
pressible Navier-Stokes equations in boundary layers is presented. The fluid domain is the 
half-space over a flat plate, and periodic conditions are applied in the horizontal direc- 
tions. The discretization is spectral. The basis functions are divergence-free and a weak 
formulation of the momentum equation is used, which eliminates the pressure term. An 
exponential mapping and Jacobi polynomials are used in the semi-infinite direction, with 
the irrotational component receiving special treatment. Issues related to the accuracy, sta- 
bility and efficiency of the method are discussed. Very fast convergence is demonstrated 
on some model problems with smooth solutions. The method has also been shown to 
accurately resolve the fine scales of transitional and turbulent boundary layers. 


1. Introduction 

The state of the art in “direct” and “large-eddy” numerical simulations of turbulent flows 
is reviewed by Rogallo &: Moin ( 1984 ). These simulations, in which the energy-containing 
structures are computed, have contributed to the understanding of several turbulent flows: 
isotropic turbulence, homogeneous turbulence with mean strain, shear or rotation, plane 
channel and curved channel flows. The results can also be used to accurately calibrate 
existing turbulence models, or suggest new models. While limited to simple geometries, 
relatively low Reynolds numbers and small statistical samples, direct numerical studies 
have the advantage of providing much more complete information than experiments can 
provide. Numerical simulations are also useful in the study of transition (Orszag & Kells 
1980 , Wray &: Hussaini 1984 ), even though the final stages of transition have not yet been 
computed for lack of numerical resolution. 

Both the numerical methods and the “supercomputers” are rapidly improving. The nu- 
merical methods used are in general spectral, at least in the periodic directions. Spectral 
methods are the most efficient in terms of resolving fine structures with a given number of 
degrees of freedom (Gottlieb &: Orszag 1977 ). They also have very high, “infinite-order” 
formal accuracy. Homogeneous turbulence was first simulated, with periodic conditions 
and Fourier series in all directions. This method is well established and major improve- 
ments are unlikely to be made. In contrast, many methods are competing for the straight- 
channel flow (Orszag & Kells 1980 , Moin Sz Kim 1982 , Kleiser 1982 , Moser et al. 1983 ). 



A few methods exist for free mixing layers (Cain et al. 1984), curved channels (Moser et 
al. 1983, Marcus 1984), pipe flow (Leonard & Wray 1982), and boundary layers (Orszag 

Patera 1983, Wray & Huss'aini 1984). 

This paper is the first of three on the direct simulation of boundary layers. Transition 
simulations are presented in Part 2, and simulations of fully developed turbulence and 
relaminarization are presented in Part 3. The purpose of Part 1 is to present a spectral 
numerical method designed specifically for the simulation of boundary layers, to discuss 
certain issues regarding accuracy and efficiency, and to report on tests the method was 
subjected to. Two major items are the weak formulation, which applies to other flows, 
and the choice of the mapping and basis functions, which is specific to the boundary layer. 
Preliminary results were presented by Spalart (1984a, 1984b) and Spalart & Leonard 
(1985). In the design of major computing projects, considerations of physical relevance 
and of computational efficiency are taken into account and sometimes conflict. For this 
project, periodic conditions were adopted in both the spanwise and streamwise directions. 
Such conditions are optimal in terms of efficiency and will be taken for granted in this 
paper. Their physical implications will be discussed in Parts 2 and 3. 


2. Weak formulation. 

Leonard (1981) presented a spectral method based on divergence-free basis functions 
and a weak formulation of the incompressible Navier-Stokes equations. He mentioned the 
following advantages over previous methods: exact treatment of the continuity condition 
and boundary conditions, simpler time-advance scheme (the pressure term is eliminated), 
and lower memory requirements. The method has most of the advantages of the vorticity- 
stream function methods, but by working with the velocity instead of the vorticity it 
avoids the well-known problem of having to prescribe vorticity boundary conditions. A 
weak formulation is used because it is consistent with the use of constrained basis functions, 
not because weak solutions are expected. The method was successfully applied to a circular 
pipe by Leonard & Wray (1982) and to straight and curved channels by Moser et al. (1983). 

The present method uses Leray’s weak formulation of the equations, which is slightly 
different from Leonard’s. Leray’s formulation was presented in 1933 and has been used 
extensively for theoretical studies of the equations (Ladyzhenskaya 1969, Heywood 1980, 
Temam 1983) and with finite-element numerical methods, but apparently not with spectral 
methods. The relative merits of these two formulations will be discussed. 

2.1. Leonard’s formulation 

The incompressible Navier-Stokes equations are written 

v.u = o, (1) 

U t + U.VU = -Vp+t/V 2 U + f, (2) 

where U = ( u,v,w ) is the velocity vector, V the gradient operator, t the time, p the 
pressure, and u the kinematic viscosity. The density, being uniform, can be set to 1 and 
omitted. The term f may be either a forcing term, or an extra term (function of U). Such 
terms will appear in the equation (see Parts 2 and 3). However they do not interfere with 
the formulation, and can be omitted for the moment. 
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The boundary conditions at a solid wall (y — 0 in this case) are the viscous conditions 

U(x,0,2) = 0. (3) 

If the domain is infinite, a “free-stream velocity” U^ is prescribed. In the present case it 
is approached for large values of y 


lim U = Uoo. (4) 

>oo 

The rate at which U tends to Uoo may have to be prescribed also; existence and uniqueness 
questions in infinite domains are delicate (Heywood 1980). 

The momentum equation (2) governs the time-evolution of the velocity field under the 
constraints of the continuity condition (1) and the boundary conditions (3,4). These 
constraints are linear and time-independent; they define a subspace of the vector space 
of all possible velocity fields. The idea is to restrict the search for a solution to that 
subspace from the onset. Numerically, the first advantage is that by working inside a 
smaller space one has fewer degrees of freedom to follow for a given level of resolution: two 
per mode, instead of four when the three velocity components and the pressure are used 
(Leonard 1981). 

Once the continuity condition (l) has been applied the description acquires a non-local 
character; the different velocity components, at different points, are coupled. This re- 
flects the instantaneous interactions, through the pressure, that characterize incompress- 
ible flows. It is then natural to apply the momentum equation, not locally and separately 
for the three velocity components, but in a global sense. This is done by taking the dot 
product of the vector equation (2) and suitably chosen vector “test” functions, and inte- 
grating over the entire domain. The inner product < Ui,U 2 > of two vector functions is 
defined by 

<U 1 ,U 2 >= 

Let V be a test function. The product of (2) with V is 

< u t , V > + < U.VU,V >= - < Vp,V > + < t/V 2 U,V > . (6) 

Leonard used test functions V that satisfy the divergence- free condition (1) and have zero 
normal velocity at the walls (they satisfy only the y component of (3)). An integration 
by parts then shows that < Vp, V >= 0; the pressure term is eliminated. Actually, any 
gradient-type term is orthogonal to V in the sense of Considering the well-known 

identity U.VU - V(|U| 2 /2) - U x w, where to is the vorticity vector, one has the choice 
of keeping U.VU or just — U x to for the transport term. Leonard arrived at the equation 

<U,V> i -<Uxw,V>=i/<V 2 U,V>. (7) 

This equation will be advanced step by step in time, and no constraints need to be applied. 
While the system (1-4) was an initial-boundary- value problem, (7) is an initial-value prob- 
lem, and is easier to integrate in time. In addition (l) and (3) are satisfied exactly, which 
is not always the case with other methods. 


/// 


Ux.U 2 dx dy dz. 


( 5 ) 
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2.2. Leray’s formulation 

In Leonard’s formulation, the test-function space is larger than the basis-function space, 
because it is defined by fewer constraints (the basis functions satisfy all three components 
of (3)). The test-function space is only “slightly” larger, in the sense that the basis-function 
space is dense in it for the norm defined by (5). However, in a numerical context, when 
finite-dimensional subspaces are used, the difference is significant. In Leray’s formulation, 
the test functions satisfy both the normal and the tangential components of (3); the two 
spaces are the same. One has a true least-squares method. Leray’s test-function space has 
the advantage of being the smallest set of test functions sufficient to govern the evolution of 
U; within that space the bilinear symmetric form < .,. > is positive-definite. Numerically, 
the advantage is that the choice of the test functions will be unique once the basis functions 
have been chosen, and that the matrix representing < .,. > will be symmetric. 

The extra constraint on V also allows one to integrate the viscous term by parts, and 
to arrive at 

<U,V> t — <Uxw,V >=-!/.< VU,VV>. (8) 

Thus the viscous operator is also symmetric, and is negative-definite. This has some 
computational advantages (lower storage, use of the Cholesky decomposition, simultane- 
ous diagonalization of the two operators), and it guarantees that the numerical Stokes 
eigenvalues will be real and negative. 

Another issue is the conservation of energy. The contribution of the transport term to 
the global kinetic energy < U,U > /2 of the numerical solution should, ideally, be 0. 
With Leray’s formulation (the expansion and test functions being the same) U can be 
constructed as a linear combination of test functions. Thus if (8) is satisfied for every V, 
one also has, by linear combination, 

-<U,U> t -<Uxw,U>=-i/<VU,VU>. (9) 

2 

The contribution of the transport term, < U x w,U >, is 0 since (U x u) and U are 
orthogonal at every point. Thus the global energy of an inviscid flow is conserved by the 
spatial discretization. Errors in the time integration may still affect the energy, so that 
one sometimes speaks of “semiconservation” of energy. However, this semiconservation 
property and the use of an adequate time-integration scheme will ensure that the numerical 
solution remains stable. Equation (9) also shows that the contribution of the viscous term 
is negative, as expected. In Leonard’s method, semiconservation of energy does not seem 
to hold; however, he does not report any stability problems. 

The disadvantage of Leray’s formulation, in a numerical context, is that Chebyshev 
polynomials cannot be used to treat the wall-bounded direction, for the reason that these 
polynomials are orthogonal with respect to a weight function, (1— y 2 ) -1 / 2 , which is singular 
at the wall. Moser et al. used Chebyshev polynomials to construct the basis functions and 
the same polynomials, multiplied by (1 — y 2 ) -1 / 2 , to construct the test functions. Thus 
the two families are orthogonal to each other with respect to the < .,. > product. The test 
functions are less regular than the basis functions near the wall. In Leray’s formulation (the 
basis functions and test functions being the same) one would have to split the Chebyshev 
weight function between the basis and test functions. As a result, the functions would have 
infinite derivatives at the wall, which must be avoided. It is essential for the accuracy of 
the spectral method to have regular basis functions (Gottlieb &; Orszag 1977). 
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This is the main reason why Jacobi, instead of Chebyshev, polynomials are used here. 
Leonard & Wray (1982) used Jacobi polynomials partly because they allowed a proper 
treatment of the polar-coordinate singularity. Numerically, Chebyshev polynomials would 
be preferable because “fast” Chebyshev transforms can be used. With these fast trans- 
forms, the number of operations the computer has to perform is of the order of N log(iV), 
instead of N 2 for a standard transform. The penalty becomes heavy when N, the number 
of polynomials, exceeds about 50. Thus Jacobi-based methods like Leonard and Wray’s 
and the present one are competitive for flows with only one wall, when the value of N is 
typically between 20 and 40, but are less attractive for channel flows, when N can be as 
high as 128. There is need for a “fast Jacobi transform.” 

Leonard’s and Leray’s formulations are elegant and lend themselves to the design of 
numerical methods with attractive properties, as will be discussed. On the other hand, the 
simplicity of the domain shape and the periodicity in two directions are used extensively. 
It remains to be seen how conveniently this formulation can be applied to more complex 
domains or to zonal methods. 


3. Numerical Method 

The discretization is spectral in space and of finite-difference type, with second-order 
accuracy, in time. In addition to their high accuracy, the global character of spectral 
methods seems to make them especially appropriate for the simulation of incompressible 
flows, since all the points in the domain are instantly coupled by the pressure interactions. 

3.1. Construction of the vector basis functions 

The three-dimensional vector basis functions are constructed as products of 
one-dimensional functions. In the x and z directions, parallel to the wall, periodic condi- 
tions are applied and Fourier series can be used. Let k = ( k x ,k z ) be the wave-vector in 
the horizontal plane and k denote its norm, \/k,' 2 + k 2 . Let U(A; X , y, k z ) denote the Fourier 
transform of U(x, y, z). Equation (l) becomes 

ik x ti + v y + ik z w = 0, (10) 

where i 2 = —1. Leonard & Wray (1982) introduced a decomposition into “+” and ” 
modes to simplify this equation. The + modes contain the horizontal component of velocity 
parallel to the wave-vector and the vertical component, and the — modes contain the 
component orthogonal to the wave- vector: u + = ( k x u + k z w)/k, u~ = ( k z ii - k x w)/k. 
Equation (10) becomes 

iku + + v y — 0. (11) 

The — modes are not involved in the continuity equation; they only have to satisfy the 
boundary conditions. They are decoupled from the + modes in both linear operators in (8). 
The + modes are two-dimensional modes in the vertical plane containing the wave-vector; 
the rotation amounts to Squire’s transformation (Schlichting 1979). 

Thanks to the Fourier transforms and to Leonard and Wray’s +/— decomposition the 
problem of defining vector basis functions is reduced to the problem of defining one- 
dimensional scalar functions with a single zero at y = 0, for u~ , and functions with a 
double zero at y — 0, for v. u + is then given by (11). The behavior of the functions as 
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y — > oo, however, deserves special attention. Leonard & Wray (1982) and Moser et al 
(1983) treated finite domains. 

3.2. Choice of the mapping 

Orthogonal polynomials make convenient “building blocks” for numerical methods. As 
shown by Gottlieb Orszag (1977) using polynomials directly in terms of y over [0, oo] is 
inefficient, essentially because all polynomials diverge as y —> oo. Grosch &: Orszag (1977) 
showed that the procedure of mapping the infinite interval [0, oo] into a finite interval and 
using polynomials in terms of the new variable was more efficient than other procedures, 
including the truncation of [0, oo] to a large but finite value of y. Two obvious candidates 
are the exponential mapping and the algebraic mapping. Orszag (1976), Grosch &; Orszag 
(1977), and Boyd (1982) all strongly criticized the exponential mapping and recommended 
algebraic mappings instead, in particular for the boundary-layer problem. Some boundary- 
layer results using an algebraic mapping are reported by Orszag &; Patera (1983). Cain 
et al. (1984) found that a hyperbolic-tangent mapping and Fourier series were efficient in 
treating the interval [— 00 , 00 ] in their study of free shear layers. 

The exponential mapping has certain advantages (simpler metric coefficient and differ- 
entiation operator, rapid expansion of the collocation points outside the boundary layer) 
and will be used here. A satisfactory solution was found to Orszag’s (1976) objection 
and the resulting method is believed to be more accurate, for our purposes, than existing 
methods. Once this mapping has been chosen a certain family of Jacobi polynomials lends 
itself well to the construction of basis functions that will keep the Stokes operators in (8) 
narrow-banded. 

Let the transformed variable p be defined by 

p == e~ y/y \ (12) 

The parameter j/o is a length scale and will be chosen to be of the order of the boundary- 
layer thickness, p varies over [0,1]. The derivatives of a function <f> satisfy 

*t = -zr*i- (13) 

yo 

Thus if <f> is a polynomial in terms of p, <j> y is also a polynomial, and of the same degree. 

The potential problem with the exponential mapping is the following (Orszag 1976). 
Far from the plate, the vorticity is assumed to tend to 0. The velocity field is then irrota- 
tional as well as divergence-free: it satisfies Laplace’s equation. For a Fourier-transformed 
variable, say v, Laplace’s equation becomes v yy — k 2 v = 0 and its solutions are linear com- 
binations of exp (ky) and exp(— ky). Only the solution exp (—Icy) satisfies the boundary 
condition (4). In terms of p , it is equal to p kyo . Now ky 0 , in practice, often takes values 
lower than 1, so that the derivatives of v with respect to p tend to 00 as p — ► 0. Clearly 
such a function cannot be well approximated by polynomials, so that slow convergence of 
the numerical method would result. A singularity was introduced by the mapping. The 
problem can be alleviated by making y 0 large in order to bring ky 0 above 1 or even higher, 
but then one is squandering collocation points away from the boundary layer itself. In ad- 
dition, any noninteger value of kyo results in unbounded derivatives at high enough order. 
With an algebraic mapping, the singularity in the mapping is weaker; the derivatives of 
the Laplace solution with respect to the transformed variable are bounded (they even tend 
to 0). Fast convergence results (Orszag 1976). 
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The problem is due to the relatively slow decay of the irrotational component of the 
solution. This leads to the idea of subtracting this irrotational component in the hope 
that the remainder, or “vortical component,” will be easier to approximate. Let us assume 
that as y — > oo the vorticity tends to 0 at least as fast as exp(— Ky) for some positive 
constant K : |cjj = 0(exp(— Ky)), and examine the behavior of the velocity. The vertical 
velocity satisfies the Poisson equation 

v yy - k 2 v — ikoj + , (14) 

where is the projection of the vorticity vector orthogonal to the (k,y) plane. Using 
(14) and (4) it can be shown that v has the form 

v == Ae~ ky + v', (15) 

where A is a constant and v' = O (exp{— Ky)) . The point is that in most flows of interest, 
K is much larger than k. Because of (11) and (15), u + can be split as v was in (15). 
Finally, the u~ component itself is 0(exp(— Ky)), because it satisfies u~ = iCj y /k, where 
uj y is the vertical component of vorticity. 

Let us now express the results in terms of p. We shall assume that the functions are 
smooth in terms of y , and only consider the effect of the singularity in the mapping. Since 
exp(-ify) = p Ky °, a function that is 0(exp{—Ky)) is of class C n (has n continuous 
derivatives) in terms of p for any integer n smaller than Kyo. Thus if Kyo is large, the 
spectral method can approximate v' with a finite but high order of convergence (Gottlieb 
& Orszag 1977). If the function is O (exp(-Ky)) for arbitrarily large values of K, infinite- 
order convergence is obtained. We conclude from (15) that, provided that the vorticity is 
0[exp( — Ky)) for a large value of K and that the error in the value of A (the “subtrac- 
tion error”) is small enough, one can approximate boundary-layer velocity fields with an 
exponential mapping and polynomials, and obtain high-order convergence. 

One can estimate the value of K in some relevant flows. In Stokes’ first problem and in 
Blasius flow, the mean vorticity is 0(exp{— Ky)) for arbitrarily large values of K (Schlicht- 
ing 1979). Thus the mean flow is approximated with infinite-order accuracy. In Stokes’ 
second problem and in the laminar sink flow, K is finite; K~ l is of the order of the 
boundary-layer thickness. One may improve the accuracy by choosing yo so that Kyo is 
an integer, without having to make yo too large. In solutions of the Orr-Sommerfeld equa- 
tion, K is equal to \Zk(Uoo — c)/2u, and is again much larger than k. The convergence to 
these solutions will be tested. 

To use the property expressed in (15), the simplest strategy is to add to the family of 
basis functions one extra function per wave- vector, a function that behaves like exp(— ky) 
as y — > oo, and to let the least-squares algorithm or the eigenvalue solver choose the value 
of A along with the other coefficients. No separate effort is made to reduce the subtraction 
error. Including an extra function adds to the complexity of the program, and could 
disrupt the structure of the matrices built around the regular basis functions. However, 
this disruption is minimized by the fact that the irrotational component, by nature, is 
orthogonal to all the other basis functions. Thus it does not even add any extra bands or 
lines in the matrices. Only the + modes require the extra function. One minor problem 
is that if kyo happens to be an integer, the extra function is not linearly independent of 
the regular basis functions, which makes the matrix singular. If kyo is close to an integer 
or large (about 4), the matrix becomes ill-conditioned. 


7 



A method that incorporates a Laplace solution as an extra basis function is likely to 
be more accurate for the following reason. With the other methods, the slowly decaying 
irrotational component is described by the “regular” basis functions. This means that some 
resolution (collocation points) must be available at distances from the wall of the order of 
fc -1 . With the present method, the extra function describes the irrotational component 
and no resolution by the regular functions is needed beyond the edge of the boundary 
layer (vortical region). In practice, which is equal to A/27T where A is the period in 
the horizontal direction, is much larger than the boundary-layer thickness. The present 
method may not yield as good an asymptotic rate of convergence (when a large number 
of polynomials are used and the error is extremely small) as was shown by Boyd (1982). 
However, with a moderate number of polynomials it is significantly more accurate. When 
solving the Orr-Sommerfeld equation, a relative error of 10~ 6 in the phase velocity is 
obtained with about 25 polynomials. Orszag (1976) needed 42 Chebyshev polynomials to 
achieve the same accuracy (better accuracy could be obtained with the algebraic mapping 
by using only the odd Chebyshev polynomials, see Spalart 1984a). 

A splitting into irrotational and vortical components may be beneficial even in a spectral 
method using an algebraic mapping, or in a finite-difference method, because it would allow 
one to bring almost all the points back into the boundary layer itself (in the finite-difference 
case, one would apply the condition v y = — kv at the edge of the domain, instead of v = 0). 
The same ideas would apply to other unbounded flows; for instance, flows around solid 
bodies or vortex rings. In summary, if a flow has a region where the velocity field is both 
incompressible and irrotational, this region can and should be treated in a very economical 
fashion. 

3.3. Choice of the polynomials 

Let us consider the family of shifted Jacobi polynomials G n of degree n, defined on [0, 1] 
by G n (l) > 0 and 

l GMG m Mridv = . (16) 

where 6 nm is the Kronecker symbol. See the appendix for other properties of the G„’s. 
We define the two families h n and g n by 

Mh) = (1 - r)hG n {rj), g n (r]) = (l - r)) 2 rjG n (rj). (17) 

The h n polynomials have a single zero at the wall and will be used for the g n polyno- 
mials have a double zero and will be used for v. Both h n and g n tend to 0 as y tends to 
+oo. As shown in the appendix, h n , g n and the required derivatives can be expressed as 
combinations of a few polynomials of the type r)G m . As a result, the matrices involved in 
(8) are banded, with five and nine bands for the h n ’s and the g n ’s, respectively. The extra 
basis function, g~i, is defined by 

g- i(h) = 3rj kyo + ( ky 0 - A)rjG 0 + (1 - ky^ijGx. (18) 

It has a double zero at the wall and is orthogonal to the regular basis functions, except 
for the few that involve Go or G\. Figure 1 is a plot of the first few g n functions versus y, 
with kyo = 0.4. The extra function jy_j stands out for large values of y. As n increases g n 
has more and more zeroes within [0, oo]; the position of the zeroes roughly indicates the 



region where good resolution is available. The collocation points are indicated by crosses 
on the axis. They are seen to cluster near the wall, and to be very sparse for y beyond 
about 5 y 0 . 

The mean component, with k = (0,0), requires special treatment. Both horizontal 
components are treated as — modes, using the h n functions, and v = 0 (Moser et al. 
1983). In addition, to obtain different values at y — 0 and for y — > oo in accordance with 
(3) and (4), the function Z7oo(l — r?) is added to the expansion in the k = (0,0) mode. This 
requires a simple correction in the viscous term. 

3.4. Accuracy tests 

The accuracy of the method was first tested by approximating typical mean velocity 
profiles: the solutions to Stokes’ first problem, the Blasius equation, and the laminar sink 
flow. The reference solution for the Blasius equation was obtained by fourth-order Runge- 
Kutta integration with a large number of steps. The relative error in the slope at y = 0 
is plotted in figure 2 as a function of the number of polynomials. Very fast convergence is 
observed for the first two cases; the convergence is only moderately fast for the sink flow, 
which confirms the analysis. These tests did not involve the extra function. 

The Orr-Sommerfeld equation was also solved; the extra function is now involved. Figure 
3 is a plot of the spectrum for the Blasius flow at its critical condition ( R$ * — 519.06381 and 
k6* = 0.303773) with 26 and 53 polynomials. The complex velocities (c r ,c t ) are plotted; 
c r gives the phase velocity and c t gives the growth rate. The principal eigenvalue is 
neutrally stable. Its convergence is also described in figure 2 and is very fast. This indicates, 
empirically, that even with the simple strategy that was adopted, the subtraction error is 
very small. The best value of yo is about 26“ . The convergence without the extra function 
y_i is also plotted for comparison and is much slower, as one would expect. 

In addition to a few discrete eigenvalues, the numerical spectra in figure 3 show a string 
of eigenvalues that starts at the point with complex velocity (1,0) and extends to large 
negative values of the imaginary part. As the number of polynomials is increased this 
string becomes denser and slowly converges to the c r — 1 axis. Grosch &: Salwen (1978) 
showed that the exact spectrum includes a continuous line on that axis, and that the 
corresponding eigenfunctions behave like sine waves as y — > oo. Such functions cannot be 
well approximated by the basis functions that were chosen, since these tend to 0 as y — » oo 
(figure 1). This explains why the convergence to the continuous part of the spectrum is so 
slow. These eigenfunctions are not relevant since they do not satisfy (4). In any case the 
representation of the continuous spectrum, while not very accurate, is not abnormal. 

3.5. Gaussian quadrature and aliasing removal 

The solution of the momentum equation (8) requires the evaluation of various integrals. 
Integrals of the type < U, V > and < VU, VV > are computed exactly using the def- 
inition of g n and h n (17), and the properties of the G n polynomials (16). The integral 
<Uxw,V > is evaluated by Gaussian quadrature in real space, which is computationally 
more efficient than a convolution sum in spectral space. In each Fourier direction the 
quadrature points are equally spaced and the number of points is 3/2 times the number of 
wave-numbers, so that the quadrature of the product (U x w).V is exact. In the polynomial 
direction the quadrature points and weights are computed using the method described by 
Golub & Welsch (1969). Again, the number of points is 3/2 as large as the number of 
polynomials and the quadratures are exact. The aliasing errors are removed. It is essential 
to use the (XJ x w) form of the transport term; the other expression, (U.VU).V, would 
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not tend to 0 fast enough as y — > oo to be integrated properly. 

The use of the “3/2 rule” is relatively expensive and it is a matter of controversy whether 
aliasing removal is actually necessary. Tests were conducted, in which a fully developed 
turbulent field was produced using the program with de-aliasing (Spalart & Leonard 1985). 
The integration in time was then continued without de-aliasing (which reduced the com- 
puter work per time step by about 40%). However, the turbulent energy quickly started 
to decay and after a few hundred time steps, the turbulence was clearly “dying.” In some 
applications of spectral methods, the solutions may be resolved well enough for aliasing 
errors to be harmless. In the present case de-aliasing is necessary. It is not clear how alias- 
ing errors disturbed the flow. One of the first symptoms was that the energy of the high 
spanwise wave-numbers increased significantly when compared with the initial spectrum. 
The problem cannot be that the global kinetic energy is affected, because one can show 
that it is semiconserved even with aliasing; it is more likely that the partition of the energy 
between low and high wave-numbers is affected. 


3.6. Time integration 

Equation (8) is integrated in time using a hybrid scheme. The viscous term 
—v < VU, VV > is better treated by an implicit scheme, in order to obtain stability 
in spite of its “stiffness.” The second-order-accurate Crank-Nicolson scheme is used. On 
the other hand, the transport term < U x w,V > is not easily treated by an implicit 
scheme, because it is hard to linearize with a spectral method. In addition, it is less stiff 
than the viscous term since it involves only one first derivative. This leads one to treat 
the transport term by an explicit scheme. A low-storage third-order-accurate scheme was 
derived by A. Wray (personal communication 1983). This scheme uses three substeps; 
the first one amounts to an Euler explicit step, while the second and the third amount to 
Adams-Bashforth steps, but with uneven substeps and weighting coefficients. The global 
scheme is equivalent to the classical third-order Runge-Kutta scheme when applied to lin- 
ear equations; the stability limit expressed as a Courant number is \/3. In practice, the 
Courant number is evaluated as 


27r , \u\ Ivl IwL 


(19) 


where Ax, Ay, and A z are the spacing between collocation points. This definition of the 
Courant number is the correct one for the Fourier directions. The definition in the y direc- 
tion was done purely by analogy. During a simulation the time step At is automatically 
adjusted so that the maximum local Courant number is \/3. Typically, because of the vis- 
cous dissipation and the intermittent character of the velocity field, numerical instability 
does not occur until the Courant number limit is set to at least 2.3. The extra term f in 
(2), if present, is treated by the explicit scheme along with the transport term. This term 
is small and does not cause any stability problems. 

Since the treatment of the viscous term is second-order accurate, the overall accuracy is 
only of second order. The third-order scheme does not improve the formal accuracy, but 
in practice it significantly improves both the accuracy and the stability (as compared with 
the commonly used Adams-Bashforth second-order scheme) and its numerical dissipation 
is moderate. An additional improvement is gained by doing the time integration in a 
Galilean reference frame that is “sliding” at a velocity c with respect to the wall. This 
reduces the value of |u| in (19), which allows longer time steps and improves the accuracy. 
For transition simulations (see Part 2), c is chosen to be the phase velocity of the primary 
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Tollmien-Schlichting wave. Thus the wave is standing in the sliding frame, and is followed 
very accurately. For turbulence simulations (see Part 3), c is taken as U^/2. There is not 
as much gain, because in these simulations A y is so small that the main contribution to 
the Courant number comes from the v component. 

3.7. Resources 

The program runs on a Cray X-MP computer and is written in the VECTOR AL lan- 
guage (A. Wray, personal communication 1983). The algorithm alternates between two 
“passes” through the data base. During the “horizontal pass,” the Fourier transforms 
and the U x w products are computed (x, z) plane by (x, z) plane. During the “vertical 
pass,” the Jacobi transforms, the Stokes matrices, and their Cholesky decomposition are 
computed ( k x ,y ) plane by (k x ,y) plane. This allows one to use secondary storage and 
access only a few planes at a time. About 40% of the computing time is used by Jacobi 
transforms, 40% by Fourier transforms, and 20% by the matrix inversions. Transition 
simulations, with 16 x 40 x 48 collocation points in x, y and z, use about 3 x 10 5 words 
of storage and 2.7 seconds of CPU time per complete time step. Each run takes a few 
hundred time steps. Typical turbulence simulations, with 192 x 45 X 128 points, use about 
10 6 words of core memory, 1.7 x 10 6 words of secondary storage (two words being packed 
into one before writing to disk), and 18 seconds per step. A fully developed turbulent field 
is obtained in a few thousand steps. 

Since the turbulence is highly inhomogeneous and anisotropic, choosing the number of 
points to be used in each of the three directions is not trivial. Given a total number 
of points, one would like the simulation to be balanced, and not to see the accuracy be 
limited by one of the directions while another one is overresolved. One may examine the 
one-dimensional spectra, but spectra represent a large amount of information and one 
would like to have a simple, objective indicator of how good the resolution is in each 
direction. In addition, spectra are not available in the polynomial direction. 

One possibility is to consider how well correlated the solution is between a grid point 
and its immediate neighbors. For small separations the correlation is directly related to 
the Taylor micro-scale A, so that it would be roughly equivalent to consider the ratio A/A. 
A is the distance between grid points. Let C x be defined by 


C x (y) 


U'(x).U'(x + Ax) 

lU'P 


( 20 ) 


where U' is the fluctuating velocity vector, and C y and C z be defined similarly. The bars 
denote averages taken in the homogeneous directions, x and z , so that C x , C y , and C z are 
functions of y. By comparing these three functions directly, one is implicitely assuming 
that the “resolving power” per grid point is the same in all directions, which is plausible 
with a fully spectral method. With a mixed spectral/finite-difference method, for instance, 
one would have to “penalize” the finite-difference direction. 

The correlations in a typical turbulent simulation are shown in figure 4. The closer 
the correlation is to 1, the better the resolution is. The values of the three correlations 
are comparable; the simulation was reasonably well balanced. As expected, C x and C z 
worsen near the wall, because the length scales decrease while Ax and A z do not. On 
the other hand, C y improves near the wall, presumably because Ay decreases faster than 
the length scales of the flow. This suggests that the Jacobi polynomials gather more than 
enough points near the wall. Ideally C y would be independent of y. This also suggests 


11 



that Chebyshev polynomials, which gather even more points near the wall, may be less 
efficient than Jacobi polynomials for the flows we are simulating. 

A comparison of the two-point correlations cannot give a rigorous proof that the resolu- 
tion is optimally distributed, but it gives a simple and convenient estimate of the accuracy, 
and helps one choose values for the major numerical parameters: N x , N y , N z , and y 0 . The 
correlations are responsive to changes in resolution; refining the grid by a factor of 2 re- 
duces the “lack of correlation” 1 - C roughly by a factor of 4. The correlation curves 
may also be used to estimate the relative accuracy of separate simulations, for instance at 
different Reynolds numbers. 


The author had useful discussions at NASA Ames Research Center with Drs. J. Kim, 
A. Leonard, P. Moin, R. Moser, R. Rogallo and A. Wray. He thanks Dr. D. Jespersen for 
reviewing the manuscript, and Mr. H. Lomax for his teaching and his support. 
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Appendix 


This appendix outlines the construction of the basis functions. The properties of the 
Jacobi polynomials that were used are the following. The two properties G n (l) > 0 and 
(16) fully define the G n family. These polynomials satisfy G n (l) = 1 and are described by 
Abramowitz & Stegun (1972) (with a slightly different normalization). The actual building 
blocks are the polynomials r)G n , which satisfy 

l = (A1) 


since dy = —y 0 dy/y. Abramowitz and Stegun also give the recursion relationship, 


yG n — 


n + 2 2 (n + l) 2 „ n 

2(2n + 3) " +1 (2n + l)(2n + 3) ” 2(2n + 1) 


(A2) 


and the differentiation relationship, 


»?(i - v) 


dG n 

dy 


n{n + 2) _ n{n + 2) n{n + 2) 

2(2n + 3) n+l (2n + l)(2n + 3) ” 2(2n + 1) n_1 


(A3) 


The functions h n and g n being defined by (17), (A2) and (A3) allow one to express h n and 
g n , as well as hn s first derivative and gn s first two derivatives, in terms of the yGm’s. 
Higher derivatives are not needed when using Leray’s formulation. For instance 


kn ~ 2(2 n 


(n + 2) 
3) 


yGn+l + 


2n 2 + 4n + 1 
(2n+ l)(2n + 3) 


yG n - 


n 


2(2n+l) 


r>G n _ 


ii 


(A4) 


dh n 

dy 


1 

Vo 


+ 2 )' 


-nG. 


J- 


(n + l) 2 


(A5) 


Note that the multiplication by (1 — y) in (17), which was done to enforce the boundary 
condition, fortuitously improves the expression of the derivative as it did for Chebyshev 
polynomials in the Moser et al. (1983) method. Equations (Al), (A4) and (A5) allow one 
to compute the matrices representing the first and third term of (8) for the - modes. Both 
matrices are pentadiagonal. The algebra is very similar for the g n functions. 
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Figure 1. Basis functions <?_ i to g 3 and collocation points. 



Number of polynomials 

Figure 2. Relative error in numerical solutions, as a function of the number of polyno- 
mials. • Stokes layer; o Blasius layer; * Sink Flow; + Orr-Sommerfeld equation, using g- 1 ; 
x Orr-Sommerfeld equation, omitting g~\. 
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Figure 3. Spectrum of Orr-Sommerfeld equation. a) with 26 polynomials; 
b) with 53 polynomials. 
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